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The  maximum  likelihood  estimates  (MLEs)  and  Cramer-Rao  bounds  (CRBs) 
for  parameters  of  harmonics  in  Gaussian  noise  have  been  well  studied.  If 
the  phase  of  the  signal  is  defined  in  the  middle  of  the  time  interval  rather 
than  at  the  beginning  (as  is  more  common),  the  information  matrix  is  ap¬ 
proximately  diagonal,  and  the  formulation  and  analysis  of  the  MLEs  and 
CRBs  are  simplified.  More  significantly,  this  simple  modification  decouples 
the  estimation  of  phase  and  frequency  and  leads  to  efficient  MLE  gradient 
descent  algorithms.  In  this  report,  these  MLE  procedures  and  CRB  analysis 
are  presented  for  the  multiple  and  coupled  harmonic  case,  as  well  as  for 
colored  noise.  A  new  criterion  on  the  required  sample  size  is  presented  to 
give  uniform  bounds  on  the  accuracy  of  diagonal  information  matrix  ap¬ 
proximation;  uniform  bounds  are  needed  to  ensure  the  effectiveness  of  the 
gradient  descent  methods.  These  methods  are  demonstrated  on  real  data 
from  battlefield  acoustic  sensors,  where  they  can  be  used  to  help  identify 
targets  of  interest,  such  as  tanks  or  trucks. 
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1.  Introduction 


The  Cramer-Rao  bounds  (CRBs)  for  estimating  the  parameters  of  harmon¬ 
ics  in  Gaussian  white  noise  were  first  presented  by  Rife  and  Boorstyn  [1,2]. 
More  recently,  the  results  have  been  presented  by  Stoica  et  al  [3]  and  Porat 

[4]  and  extended  to  include  multiplicative  noise  by  Zhou  and  Giannakis 

[5] .  The  estimation  algorithms  presented  by  Stoica  et  al  [3]  do  not  exploit 
a  key  observation  made  originally  by  Rife  and  Boorstyn  [1],  and  more  re¬ 
cently  by  Porat  [4],  on  the  effect  of  defining  the  phase  of  the  signal  in  the 
middle  of  the  time  interval  rather  than  at  the  beginning,  which  is  more  com¬ 
monly  done.  Perhaps  some  reluctance  in  using  the  less  common  definition 
is  due  to  the  misperception  (as  in  Porat  [4])  that  odd  and  even  numbers 
of  samples  need  to  be  treated  separately.  With  the  phase  in  the  middle  of 
the  time  interval,  the  information  matrix  is  approximately  diagonal  as  com¬ 
pared  to  the  results  in  Stoica  et  al  [3]  and  elsewhere,  where  the  phase  and 
frequency  estimates  have  significant  correlation,  and  the  formulation  and 
analysis  of  the  maximum  likelihood  estimates  (MLEs)  and  the  CRBs  are 
simplified.  More  significantly,  this  simple  modification  decouples  the  esti¬ 
mations  of  phase  and  frequency  and  leads  to  more  efficient  MLE  gradient 
descent  algorithms.  To  ensure  the  effectiveness  of  the  gradient  descent  al¬ 
gorithms,  I  suggest  a  new  criterion  for  the  required  sample  size  that  gives 
uniform  bounds  on  the  accuracy  of  diagonal  information  matrix  approxi¬ 
mation.  In  this  report,  these  MLE  methods  and  CRB  analysis  are  presented 
for  the  multiple  and  coupled  harmonic  cases,  as  well  as  for  colored  noise, 
with  application  to  real  data  from  battlefield  acoustic  sensors.  These  ap¬ 
proaches  and  analysis  provide  valuable  insights  into  the  problem  of  target 
identification  and  the  design  of  computationally  efficient  battlefield  acous¬ 
tic  systems. 
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2.  Harmonic  Signal  Models 


The  classic  real  single  harmonic  signal  model  for  the  signal  s(t)  at  continu¬ 
ous  time  t  is 

s(t)  =  Acos(ut  +  4>)  ,  (1) 

where  A  >  0  is  the  amplitude,  u  is  the  angular  frequency,  and  </>  is  the 
phase.  The  frequency  /  =  ui/(2n)  is  often  the  natural  parameter,  but  using 
ui  helps  simplify  the  notation,  and  results  can  be  readily  applied  to  both.  It 
is  sometimes  preferable  to  write  equation  (1)  as 

s(t)  =  ucos(u>t)  +  vsin (cat)  ,  (2) 

where  u  =  Acos(4>)  and  v  =  -A  sin(<j>)  are  the  in-phase  and  quadrature- 
phase  components  of  the  signal.  The  parameters  in  equation  (2)  can  be 
transformed  to  those  in  equation  (1)  via 

A  =  (u2  +  v2)1/2  , 

4>  =  arctan(— v,  u )  ,  (3) 

where  arctan  is  the  four-quadrant  arc  tangent  function  (MATLAB  function 
ATAN2).  The  parameters  of  this  model  will  be  denoted  by  0  =  (A,  </>,  ui)  or 
6  =  ( u,v,f ). 

The  multiple  harmonic  model  generalizes  equation  (1)  as 

m 

s(t )  =  Y^Ak  cos  (ujkt  +  <j>k)  (4) 

k= 1 

and  equation  (2)  as 

m 

s(t)  =  ^2  uk  cos(ukt)  +  vk  sin (ukt)  ,  (5) 

k= 1 

so  that  the  multiple  harmonic  parameters  are  6  =  (A,  (f> ,  w)  in  equation  (4) 
and  9  =  (u,  v,  ut)  in  equation  (5).  An  important  model  that  is  a  special  case 
of  equation  (4)  is  the  coupled  harmonic  signal  model 

m 

Ak  cos(hkUt  +  fa)  ,  (6) 

fc= i 
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where  u  is  called  the  fundamental  frequency,  hk  is  an  integer  indicating  the 
/rth  harmonic  number  (assumed  known),  and  .-1  /,■  and  <pk  are  the  amplitude 
and  phase  of  the  kth  harmonic.  It  is  clear  that  equation  (6)  is  equivalent 
to  equation  (4)  with  =  hkuj,  but  the  number  of  parameters  has  been 
reduced  (2m  +  1  instead  of  3m)  and  will  be  denoted  by  6  =  (A.  <fi,  u)  or 
6  =  (u.  v,ca). 

The  data  sequence  A"  [m]  received  at  the  acoustic  sensor  (for  the  application 
presented  in  sect.  7)  is  an  equally  spaced  sampled  version  of  the  continuous¬ 
time  model: 

X[m]  =  s(t(m))  +  Z[m\  ,  (7) 

where  t(m)  is  the  time  of  the  mth  sample  and  Z[m]  is  the  cumulative  addi¬ 
tive  noise.  Without  loss  of  generality,  it  is  assumed  that  the  sampling  rate 
is  T  =  1  so  that  difference  t(m)  —  t(m  —  1)  =  1  for  all  m,  and  the  frequencies 
are  normalized  so  that  Nyquist  is  n.  In  order  to  ensure  that  the  parameters 
are  identifiable,  it  is  required  that  A  >  0  and  0  <  oj  <  tt  (0  <f  <l/2). 

Let  y  =  (yi,  ij2, . . . .  yn)T  be  a  vector  denoting  a  given  set  of  n  samples  from 
the  data  over  a  particular  time  window.  The  time  reference  in  equation  (1) 
can  always  be  redefined,  so  without  loss  of  generality,  we  can  assume  that 
yi  =  X  [<(*)].  The  model  in  equation  (7)  can  be  written  as 

x  =  s (0)  +  z  ,  (8) 


where 

Sj(ff)  =  Acos(cjt(i)  +  4>)  .  (9) 

First  consider  the  case  where  Z  is  assumed  white,  so  that  the  noise  vec¬ 
tor  z  is  multivariate  Gaussian  Nn( 0,  a2 1);  that  is,  zi,  Z2,  ■  •  ■ ,  zn  are  indepen¬ 
dently  and  identically  distributed  (i.i.d.)  Gaussian  with  variance  a2.  To  sim¬ 
plify  the  analysis,  I  assume  the  variance  is  known,  although  all  the  results 
presented  are  equally  valid  in  the  unknown  case.  In  this  case,  the  MLE 
maximizes 

log(L(y,0))  =  —  2^2  Hy  —  s(0)||2  —  n/2  log(27r<72)  ,  (10) 

which  is  equivalent  to  minimizing  the  residual  sum  of  squares 

SS(0)  —  ||y  —  s(0)|2  =  ||e||2  ,  (11) 


where  e  =  e(9)  is  the  residual  error. 

The  information  matrix  for  estimating  the  parameters  of  s(0)  is  given  by 


irm  = 

^  a2dd  89 


(12) 
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where 


(13) 


ds  1  _  <9sj 

d6_  ij  ddj 

is  a  matrix  of  first  derivatives. 

Calculations  of  the  CRBs  for  the  harmonic  signal  models  are  somewhat 
complicated  but  involve  the  straightforward  use  of  trigonometric  identities 
and  summation  formulas.  While  exact  information  matrices  can  be  readily 
calculated,  interpretation  is  only  practical  if  approximations  are  used  for 
large  sample  sizes  n.  Approximate  results  are  obtained  via  the  standard 
summation  formula  (and  its  first  and  second  derivatives),  which  for  arbi¬ 
trary  time  indices  is 

—  ^2  cos (u>t(i)  +  </>)  =  cos(<p  +  (t(l)  +  (n  —  l)/2)u;)e(u;/2)  ,  (14) 

77-  .  1 
i=l 


where 


sin(nu;) 

nsin(u;) 


(15) 


is  the  normalized  Dirichlet  kernel.  A  classic  result,  first  presented  by  Rife  and 
Boorstyn  [1]  (for  the  single-tone  complex  case)  and  also  discussed  in  detail 
by  Kay  [6],  is  the  CRB  for  the  single  harmonic  model,  where  the  information 
matrix  is  approximately 


I  (A,<t>,U>)  S3  p 


o  n  £?=1i(i)  > 

0  Etit(0  Et iKi)2. 


(16) 


where  p  =  A2 /{2a2)  is  the  signal-to-noise  ratio  (SNR).  For  Fourier  analy¬ 
sis  and  other  applications,  the  time  samples  are  taken  to  be  t(i)  =  i  —  1 
or  t(i)  =  i,  so  that  the  phase  is  associated  with  the  beginning  of  the  time 
window.  For  t(i)  =  i  —  1,  inverting  equation  (12)  leads  to  the  approximate 
CRBs 


CRB  [A]  S3 
CRB  [4>]  « 


A2  _  2 a2 
np  n  ’ 
2(2n  —  1) 
(n(n  +  1  )p  ’ 


for  phase  and  amplitude  and 


CRB[o>]  S3 


12  12 
n(n2  —  l)p  n3p 


(17) 


(18) 
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(19) 


for  frequency.  It  also  leads  to 


CRB[/]  = 


1 

(2^ 


CRB[tj]  ~ 


(27r)2n(n2  —  1  )p 


for  the  alternative  definition  of  frequency. 

As  equation  (12)  shows,  the  approximate  information  matrix  can  be  made 
diagonal  if  the  sum  of  the  time  samples  is  0.  This  can  be  accomplished  with 
t(i)  =  i  —  (n  +  l)/2  and  corresponds  to  associating  the  phase  with  the 
middle  of  the  time  window.  In  this  case. 


I  (40,  uj)  &  np 


1 

A1 

0 

0 


0  0 

1  0 

0  liz(n2  ~~  ■*•) 


(20) 


which  can  be  readily  inverted  to  give  approximate  CRBs  that  agree  with 
equation  (17),  except,  notably,  for 

CRB[<jf>]  «  —  ,  (21) 

np 

which,  as  first  noted  by  Rife  and  Boorstyn  [1]  is  the  minimum  value  for  any 
location  of  phase  and  is  equal  to  the  CRB  when  frequency  is  known. 

The  bound  for  estimating  the  phase  in  the  middle  of  the  time  window  is  a 
factor  of  2(2n  —  l)/(n+ 1),  or  roughly  four  times  smaller  than  the  bound  for 
estimating  the  phase  at  the  beginning  of  the  time  window.  To  see  why,  note 
that  the  estimate  of  phase  at  the  first  sample  and  the  phase  at  the  middle  of 
the  samples,  denoted  by  0O  and  4>l/ 2,  respectively,  are  related  by 


0o  =  0i/2  ~(n-  l)u>/2  , 

so  that  for  minimal  variance  unbiased  estimators, 

12 


Var[0o]  =  Var[01/2]  +  Var[(n  -  l)w/2]  =  ^  +  ~  4  ~  n^2 


_ _  2(2n  —  1) 

1  )p  n(n  +  1  )p  ’ 


(22) 


(23) 


which  agrees  with  equation  (17). 

The  reduction  in  variance  of  the  estimate  of  0 1/2  should  not  be  miscon¬ 
strued  as  somehow  a  better  overall  estimate  of  the  signal.  The  equality  in 
equation  (22)  (or  similar  relationships)  can  be  used  in  reverse,  to  go  from 
equation  (21)  to  equation  (17)  (or  other  CRBs),  so  if  one  is  interested  in  0o 
or  any  other  location  of  phase,  the  variance  will  increase  and  nothing  will 
change  that.  In  terms  of  ease  of  analysis  for  the  single  harmonic  case,  the 
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tradeoff  in  using  the  phase  in  the  middle  of  the  window  is  between  invert¬ 
ing  a  2  x  2  matrix  in  equation  (16)  and  evaluating  J2 1 (*)2  for  nonstandard 
time  indices  (i.e.,  t(i)  =  i-(n  + 1)/2  instead  of  t(i)  =  i  -  1).  The  preference, 
debatable  at  best,  is  certainly  not  a  big  issue.  However,  with  more  com¬ 
plicated  analysis  such  as  that  presented  by  Stoica  et  al  [3]  and  Zhou  and 
Giannakis  [5],  the  small  increment  of  simplicity  in  this  approach  is  prob¬ 
ably  welcome.  Unless  noted  otherwise,  the  results  presented  here  reflect 
the  phase  defined  in  the  middle  of  the  time  window.  Besides  easing  the 
analysis,  it  is  shown  in  section  5  that  the  diagonality  of  the  information 
matrices  leads  to  more  efficient  estimation  procedures,  including  instances 
in  which  the  decoupling  of  the  frequency  and  phase  estimates  is  crucial. 

An  approximate  CRB  for  harmonic  signal  models  can  be  made  more  ex¬ 
plicit  through  use  of  an  asymptotic  Cramer-Rao  bound  (ACRB),  which  can 
be  obtained  from  an  asymptotic  information  matrix  of  normalized  param¬ 
eters.  For  example,  in  the  case  of  a  single  harmonic,  the  ACRBs  are  given 
by 


ACRB  [A]  =  2a2  , 

ACRBM  =  -  , 

P 

ACRB  [tv]  -  —  , 
P 


and  are  obtained  by  inverting 


loo  (A,  </>,&)  =  lim  I„(n  1/2  A,  n  1/2<t>,n  3/V)  =  p 


1 

A7 

0 


1 

12 


(24) 


(25) 


where  1^  is  the  asymptotic  information  matrix.  In  this  report,  each  CRB 
given  as  an  approximation  is  with  the  understanding  of  this  relationship  to 
the  ACRB. 


Stoica  et  al  [3]  have  shown  that  all  the  cross-terms  in  the  asymptotic  infor¬ 
mation  matrix  for  the  multiple  harmonic  signal  model  involving  different 
harmonics  are  zero,  so  the  generalization  of  equation  (20)  is  still  approxi¬ 
mately  diagonal. 


I(  A,cf>,u>) 


10  0 
0  D  0 
0  0 


(26) 
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where 


D  =  diag(Af, . . . ,  A^)  (27) 

is  a  diagonal  matrix.  So  the  CRBs  for  the  multiple  case  are  approximately 
all  the  same  as  the  single  case;  that  is 

CRB[i,] 

CRBfe] 

CRBy,.] 

where  pu  =  Af./ (2a2)  is  the  SNR  of  the  fcth  harmonic. 

For  the  coupled  harmonic  case,  the  cross-terms  are  all  asymptotically  zero, 
and  the  information  matrix  becomes 


~<T 


1 


npk 

12 

n5Pk 


(28) 


/(A.  4>,u) 


n 

2a2 


I  0 
0  D 
0  0 


0 

0 


hn 2 


EfcLi  h\A2 


(29) 


so  that  the  CRBs  for  amplitude  and  phase  are  the  same  as  in  equation  (28). 
The  CRB  for  fundamental  frequency  is 

to  /  ™  \ -1 

CRB[£]«-  ,  (30) 

n  \fc= l  / 


which  agrees  with  the  result  inSwami  and  Ghogho  [7],  which  was  obtained 
by  considering  the  coupled  harmonic  model  as  a  constrained  version  of  the 
multiple  case. 
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3.  Summation  Formulas  for  Well-Resolved  Frequencies 


All  the  entries  of  the  information  matrices  for  harmonic  signal  models  in¬ 
volve  summations  of  the  form 

n 

6p(ui,u>2,<t>i,<h)  =  £  t(i)p  cos(u>ii(i)  +  di)  cos(o>2i(i)  +  ^2)  (31) 

i=l 


and 

ep{oj,  <j>)  =  J2  *(*)p  cos(ut(i)  +  <t>)  (32) 

n  i=l 

for  p  =  0, 1,  2,  which  are  related  by 

5p(u>i,LU2,4>l,<p2 )  =  ^nP+1  (ep(u  1  “  ^2,^1  -  ^2)  +  Cp(^l  +  ^2,</>l  +  ^>2))  (33) 

with  the  use  of  a  standard  trigonometric  identity.  Using  the  above  notation, 
one  finds  that 


dsT  ds 
dA{  dAj 
dsT  ds 
d<pi  dtfij 
dsT  ds 
du)i  duj 
dsT  ds 
dAi  d(f)j 
dsT  ds 
dA{  duij 
dsT  ds 
dcp-i  duJ3 


=  60(uii,U]j,di,dj)  , 

—  AiAjSo^UJi)  i Oj)  di  +  7T / 2,  (f)j  -I-  7t/2)  , 

—  AiAjS2^^i^  bJj ,  di  -I-  tt/2,  dj  vr/2)  , 

—  Aj Sq  (cjj ,  ;  di  1  dj  *1”  ^”/2)  ; 

=  Aj5i(u)i,u>j,di,dj  +  tt/2)  > 

— -  AiAj$\(uJi) ujj )  di  ^r/2?  dj  "h  ^"/^)  ? 


(34) 
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are  the  exact  terms  of  the  information  matrix  for  multiple  harmonics. 


Sq  ( hiiv ,  hj uj.  (pi ,  <pj )  , 

AjAj 5o{hiU,  hjuj ,  (pi  +  7r/2,  <pj  +  n/2)  , 

777  777 

EE  hihjAiAj52(hiUJ ,  hju> ,  (pi  +  7r/2,  <pj  +  tt/2)  , 
i=i  j=i 

^4j<5o(/izu;,  hju <pj  +  7r/2)  ,  (35) 

777 

y;  hjoo,  (pi ,  <pj  +  7r/2)  , 

777 

hjAi  ^  ^  hjAjd\  {hiLJ,  hjU),  4>i  A  tt /2,  <pj  +  7r /2)  , 

7=1 

are  the  exact  terms  of  the  information  matrix  for  coupled  harmonics. 
Approximations  of  the  information  matrices  of  harmonic  signal  models 
(co»j  =  hioj  for  coupled  harmonics)  require  that 

(i/)  (klj  ?  UJj  ;  fii  )  <f>]  )  ^  0 , 

5p(ui,Ui,<t)h<t)i)  «  inp+1ep(0,0)  ,  (36) 

for  1  <  i  <  j  <  m.  The  approximations  in  equation  (36)  are  valid  if  the 
summations  in  equation  (32)  are  near  0  for  all  combinations  of  their  sums 
and  differences,  uj,  ±u>j,  except  for  u>i  =  Uj .  If  this  requirement  is  met  (under 
some  specified  criteria),  the  frequencies  will  be  called  well  resolved.  For  a 
given  set  of  frequencies,  a  natural  question  is  how  large  n  needs  to  be  in 
order  for  them  to  be  well  resolved. 

Using  two-term  expansion  of  the  exact  information  matrix  about  the  ap¬ 
proximation  in  equation  (26),  Porat  [4]  shows  that  for  multiple  harmonics, 
the  first-order  accuracy  of  the  approximate  CRBs  in  equation  (28)  is  depen¬ 
dent  on  the  value  of  the  normalized  Dirichlet  kernel  in  equation  (14).  The 
normalized  Dirichlet  kernel  attains  a  maximum  of  1  at  w  =  0  (and  multiples 
of  7 r)  and  has  main  lobes  between  ±7r/n  and  n  ±  n/n.  A  reasonable  require¬ 
ment  for  n  is  that  these  critical  values  of  u  lie  outside  these  main  lobes;  this 
leads  to  Porat's  requirement  [4], 


Similarly, 

dsT  ds 
dAi  dAj 
dsT  ds 
d<t>i  dcpj 
dsT  ds 
du>  dw 

dsT  ds 
dAi  d(f>j 
dsT  ds 
dAi  du 

dsT  ds 
d(j>i  dw 


n  >  max 

1  <7  <7  <777 


7 r 

7 

LJi 


7 r 

7T  -  U )i 


27 r 


2ty 


|c <7  2  Cc j  27T  UJ i  UJ  j 


(37) 
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which  is  a  simple  rule  of  thumb  for  estimating  a  sufficient  sample  size. 
While  equation  (37)  controls  the  bounds  on  the  variances  of  the  param¬ 
eters  via  equation  (49),  it  is  also  important  that  covariance  terms  be  small, 
such  as  when  using  the  method  of  scoring.  The  correlations  between 
frequency  and  both  phase  and  amplitude  are  on  the  order  of  ei  and 
ei(7r/n)  «  -I/7 r  =  -0.3183  for  the  minimal  sample  size  to  satisfy  equation 
(37).  This  may  not  be  small  enough  for  some  applications,  and  equation  (37) 
would  underestimate  the  required  sample  size  in  this  instance.  I  propose  a 
supplemental  criterion  for  calling  frequencies  well  resolved  that  ensures 
that  the  summations  in  equation  (32)  are  uniformly  small. 

Equation  (14)  shows  that,  more  specifically,  the  denominator  of  equation  (15) 
must  be  large  in  order  for  e0  to  be  small.  This  motivates  calling  n  sin(w)  the 
frequency  resolution  and  defining  its  inverse. 


A(co)  = 


_ 1_ 

n  sin(w)  ’ 


(38) 


which  could  be  called  the  frequency  coarseness.  With  the  phase  defined  in  the 
middle  of  the  time  window,  t(  1)  =  -(n  -  l)/2,  and  equation  (14)  simplifies 
to 

eo(u,  <f>)  =  cos(</>)e(o;/2)  =  cos(^)  sin(nu>/2)A(u;/2)  ,  (39) 

which  is  valid  for  n,  both  odd  and  even.  A  particular  case  of  equation  (39) 
is 

n 

eo(ce,  —  7r/2)  =  ^sin(u;f(i))  =  0  (40) 

i=  1 

for  all  cj,  which  is  also  evident  by  the  antisymmetry  of  the  time  samples  and 
the  symmetry  of  the  sine  function  around  0.  Differentiating  eo(w,  <j>  —  n/2) 
leads  to 

^  sin(0)A(w/2)  (—  cos(nu/2)  +  sin(na;/2)  cos(w/2)A(a;/2))  ,  (41) 


which,  by  the  use  of  l'Hopital's  rule,  can  be  shown  to  be  continuous  at 
u>  -  0,  where  ei(0,  <fi)  =  0.  Similarly,  differentiating  e\ (u,  <j)  —  7t/2)  leads  to 


e2(w,  f)  =  i  cos(^)A(u>/2)  (-  sin(no;/2)  +  2cos(nuj/2)  cos(o;/2)A(a;/2) 

—  sin(nw/2)(l  +  cos2(nw/2))A2(ce/2)^  ,  (42) 

which  is  also  continuous  atw  =  0,  where,  for  example,  e2(0,0)  =  (1  -  n_2)/12. 
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It  is  easy  to  find  the  following  bounds, 


ko(w,0)|  <  A(w/2)  , 

\eiM) |  <  (A(ta/2)  +  A2(u/2))/2  ,  (43) 

|«2(w,^)|  <  (A(w/2)  +  2A2(ca/2)  +  2A3(w/2))/4  , 

which,  for  all  pertinent  w  >  0  and  n  already  satisfying  equation  (37),  leads 
to  the  simple  uniform  bound  on  equation  (32), 

lep(w)^)l  £  A(cj/2),  (44) 


and,  ultimately, 

|^p(wi,^2,  ^i,^2)|  <  np+1A(u;/2)  ,  (45) 

for  u>i  7^  W2  in  equation  (33).  So  perhaps  better  criteria  for  selecting  n  are  to 
require  that  A  be  small  and  to  ensure  the  uniform  bounds  in  equations  (44) 
and  (45),  which  lead  to 


A(u>)  =  min  [A(u;*),  A  flu;*  -  coj  1/2),  A((c^  +  Uj)/ 2)]  <  e 

for  some  specified  tolerance  e  >  0.  These  criteria  are  equivalent  to 

1 _ 1 _ 1 

sin(wj) !  sin(|cai  —  Uj\/2)  ’  sin((u>i  +  Uj)/ 2)  ’ 

which  is  nearly  the  same  as 

1  /  1  1  2  2 
n  >  max  -  — , - ,  - - r,  - - 

1  <i<j<m  e  \u>i  7T  —  U>i  \U)i  —  U)j\  Zn  —  Ui  —  Ulj 


1 

n  >  max  - 

l<i<j<?77,  € 


(46) 


(47) 


(48) 


when  one  uses  approximations  of  the  sine  function  near  0  and  n.  Note  that 
equation  (48)  is  a  generalization  of  equation  (37),  with  e  —  l/i r.  So  in  or¬ 
der  to  get  stricter  requirements  on  the  sample  size,  the  tolerance  should  be 
selected  so  that  e  <  l/7r. 

In  terms  of  A,  the  expressions  given  by  Porat  [4]  for  multiple  harmonics 
become 

A 

CRB[ijt]  =  -^-{l  -cos(2^)sin(nu;fc)A(a;fe)-|-o[A(u;)])  , 

CRB[<4]  =  —  { 1  +  cos(20fc)  sin(nwfc) A(cafc)  +  o[A(a>)] }  ,  (49) 

npk  1  ; 

for  amplitudes  and  phases,  and 

CRB[u>fr]  =  4?-{i  +  3cos(2<£A.)  sin(ruvk)A(uk)  +  o[A(m)]|  (50) 

»Pk  1  J 
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for  frequency,  where  o[A(u>)]  is  a  term  that  goes  to  0  faster  than  A(w). 
The  analogous  expression  for  the  CRB  between  frequency  and  phase  of  the 
same  harmonic,  for  example,  is 

CRB[<j)k,u>k]  =  cos{nuk)A(LOk)  +  o(A{u>))  ,  (51) 

nzpk 

which  suggests  the  appropriateness  of  equation  (46)  if  this  term  needs  to 
be  near  0. 

The  requirements  for  well-resolved  frequencies  generally  apply  to  the 
coupled  harmonic  signal  model,  and  expressions  like  equation  (49)  are 
readily  available  for  amplitudes  and  phases.  The  analogous  expression  of 
equation  (50)  for  fundamental  frequency  is  not  especially  illuminating,  but 
the  ratio  of  the  second  term  to  the  first  term  is  seen  to  be  on  the  order  of 
A(u>)  when  one  notes  that 


asTds  -5>A(0,0)f >M!  <  ,'AMEEW 


duu  duj  2 


i=l 


i— 1  j= 1 


=  n3 A(cj)  ^  hjA 

m 

<  mns A(u>)  ^  hf  A 


(52) 


2 

i  5 


z=l 


where  the  last  step  follows  from  the  Cauchy-Schwarz  inequality.  So  the 
first-order  approximations  for  the  coupled  harmonic  case  are  still  on  the  or¬ 
der  of  A(o?),  and  the  choice  of  n  by  equation  (46)  or  (47)  is  still  appropriate. 


4.  Maximum  Likelihood  Estimation  for  Multiple  and  Coupled 
Harmonics 


Since  equations  (6)  and  (1)  are  special  cases  of  equation  (4),  we  present  the 
MLE  procedures  in  terms  of  the  multiple  harmonic  case.  The  gradient  de¬ 
scent  methods  for  equation  (6)  are  sufficiently  distinct  and  will  be  treated 
separately.  With  equation  (5),  the  multiple  harmonical  parameter  estima¬ 
tion  problem  can  be  written  as  a  linear  model 

y  =  X/3  +  z  ,  (53) 

where  X  =  X(u>)  is  an  n  x  2m  matrix  and  the  parameter 


P  = 


u 

V 


(54) 


is  a  2m  x  1  vector.  One  can  write  X  =  [C  S],  where  C  and  S  are  n  x  m 
matrices  given  by 


C  ij  =  cos[uq-f(i)] 

Sj j  =  sin[wjf(i)]  ,  (55) 

so  that  equation  (53)  becomes 

y  =  Cu  +  Sv  +  z  (56) 


in  terms  of  u  and  v.  For  Gaussian  white  noise,  the  MLE  of  (3  for  fixed  fre¬ 
quencies  to  is  also  the  least-squares  estimator  (LSE),  given  by 

P  =  PH  =  (XTX)“1XTy  (57) 


when  the  standard  method  of  least  squares  is  used.  The  MLEs  of  u  and  v 
are  readily  found  with  equation  (54),  and 


A-k  —  (u2k  +  vl)1^2  , 

4>k  =  arctan(— %,  uk) 


(58) 


are  the  MLEs  of  the  amplitude  and  phase.  Finding  the  MLE  of  9  requires 
some  sort  of  search,  either  grid  or  gradient  descent,  to  find  the  minimum 


13 


value  of  equation  (11).  So,  for  moderate  m,  finding  the  MLE  of  the  har- 
monical  parameters  is  potentially  computationally  expensive  because  of 
the  matrix  inversion  in  equation  (89).  Fortunately,  with  the  phase  defined 
in  the  middle  of  the  time  window,  XTX  is  block  diagonal  and  it  decouples 
the  estimation  of  u  and  v.  Moreover,  XrX  is  approximately  diagonal,  and 
accurate  estimates  can  be  obtained  without  inversion. 

In  terms  of  C  and  S, 


xTx  =  [  cT°  sT°  1  =  [  °Tc  ° 

CTS  STS  0  srs 


(59) 


The  estimation  of  u  and  v  decouple  into 

u  =  (CTC)~1CTy  , 

v  =  (STS)_1STy  ,  (60) 


which  is  computationally  more  efficient  than  using  equation  (57)  if  the 
phase  is  not  defined  in  the  middle  of  the  time  window  (inverting  two  m  x  m 
matrices  is  better  than  inverting  one  2  m  x  2m  matrix).  For  single  harmonics, 
the  exact  MLE  can  be  easily  found,  because 

CTC  =  |(l  +  eM), 

STS  =  £(l-eM),  (61) 

where  e(u>)  is  given  by  equation  (15).  So 


u  = 

V  = 


1  z  r^T 

TT7RnC  y' 
1  2 

1 - TT_S  y 

1  —  e(cj)  n 


(62) 


is  the  exact  MLE  from  equation  (60). 

For  multiple  harmonics,  the  exact  entries  of  the  matrices  to  be  inverted  are 


[cTc]^  =  (5o(cai,u)j,0,0)  , 

STs]  =  So(u>i,u)j,-Tr/2,—n-/2),  (63) 

which  for  well-resolved  frequencies  leads  to  the  well-known  result 

XrX  »  |l,  (64) 
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with  the  use  of  the  approximations  in  equation  (36).  So  accurate  approxi¬ 
mate  MLE  solutions  without  matrix  inversion  for  u  and  v  are 

P  ~  -XTy  ,  (65) 

n 


which  in  the  single  harmonic  case  is  equivalent  to  setting  e(w)  to  0  in  equa¬ 
tion  (62).  An  alternative  approximate  MLE  for  the  multiple  case  is  to  use 
equation  (62)  for  each  estimate  separately. 


uk 


h 


1  2 
1  A  c(ca)  vt 


_L _ 

1  -  e(w) 


(66) 


where  C  =  [Cls  C2, . . . ,  Cm]  and  S  =  [Si,  S2,  •  •  • ,  Sm].  The  estimates  in 
equation  (66)  correspond  to  setting  all  the  off-diagonal  terms  in  CTC  and 
STS  equal  to  zero  in  equation  (60). 

The  MLE  of  u  and  v  can  also  be  efficiently  and  simply  estimated  with  the 
fast  Fourier  transform  (FFT)  of  y.  The  continuous  Fourier  transform  (CFT) 
of  y,  denoted  by  Y (/),  is 

nf)  =  jre-2*mi)yi  (67) 

i= 1 


for  0  <  /  <  1/2.  Usually  t(i)  =  i  —  1  in  equation  (67),  but  for  the  com¬ 
parison  of  different  phase  estimates,  it  is  important  to  stay  consistent.  In 
practice,  one  can  easily  add  a  phase  correction  to  equation  (67)  to  account 
for  any  differences  in  the  time  indexes.  The  FFT — or  more  precisely,  the  dis¬ 
crete  Fourier  transform  (DFT) — calculates  equation  (67)  only  at  the  center 
frequencies  Fk  =  k/n  of  the  Fourier  bins  for  k  =  0, 1, . . . ,  (n  -  1).  For  each 
u>  =  2ni,  one  can  obtain  approximate  estimates  using  the  FFT,  by  rounding 
the  frequencies  to  the  nearest  frequency  bins: 

uk  =  ^Re{F[round(n/fc)/n]}  . 

vk  =  -^Im|y[round(n/fc)/n]}  ,  (68) 

where  round  (x)  is  the  nearest  integer  to  x.  This  procedure  is  equivalent 
to  equation  (57),  with  the  harmonic  frequencies  rounded  off  to  the  nearest 
Fourier  bin,  which  makes  X 'r X  =  (n/ 2)1  (exactly).  The  MLE  of  the  ampli¬ 
tude  (of  the  harmonic  at  the  nonrounded  frequencies)  can  be  obtained  with 
equation  (58),  but  to  account  for  using  the  center  frequency  of  the  Fourier 
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bins  rather  than  the  exact  frequencies,  a  better  estimate  of  the  phase  (i.e., 
less  biased)  is 


<t>k  =  arctan(— vk,  uk )  -  nfk  -  round (nfk) 


n 


7 r- 


n 


(69) 


which  can  be  a  correction  of  almost  tt/2  for  harmonic  frequencies  that 
are  half  a  bin  away  from  the  center.  Note  that  the  approximate  estimate 
equation  (65)  is  equivalent  to  use  of  the  CFT  in  equation  (67)  rather  than 
the  DFT,  which  leads  to 


Uk 


h 


-Re 

n 

--Im 

n 


Y{fk) 
Y(fk ) 


(70) 


this  is  simply  equation  (68)  without  rounding. 


5.  Gradient  Descent  Methods  for  Finding  MLE 


The  MLE  of  the  harmonic  parameters  can  be  found  by  gradient  descent 
methods.  Since  for  fixed  frequencies,  the  MLE  of  the  other  parameters  can 
be  obtained  with  equation  (60),  these  gradient  descent  methods  are  needed 
to  find  the  MLE  of  the  frequencies  cu  or  the  fundamental  frequency  u.  Gra¬ 
dient  descent  methods  can  be  implemented  directly  in  terms  of  the  ampli¬ 
tudes  and  phases,  but  this  does  not  take  advantage  of  the  linearity  of  the 
parameterization  in  equation  (5).  The  Newton-Raphson  method  uses  an 
initial  estimate  Go  and  iteratively  updates  the  estimate  via 


d2log[L(y,fl)] 

dOdeT 


dlog[T(y,<9)] 

de 


which  requires  the  inversion  of  the  Hessian  of  the  log-likelihood  function. 
The  gradient  vector  and  the  Hessian  matrix  in  equation  (71)  can  be  written 

as  _ 

01og[L(y,0)]  _  1  dsT^ 


8GdO J 

where  the  second  component  is 


■>(»)  =  is 


^d*Si 


o  ^z 

dd2 


and  I  is  the  information  matrix  defined  by  equation  (12).  An  alternative  to 
the  Newton-Raphson  method,  called  the  method  of  scoring,  is  to  replace  the 
Hessian  by  its  expected  value;  this  method  of  scoring  is  often  preferable 
because  it  involves  less  computation  and  has  stabilizing  properties.  From 
equation  (74),  £[J(0)]  =  0,  and  the  expected  value  of  the  Hessian  is 

\  dGdGT  j 

so  that  replacing  the  Hessian  in  equation  (71)  by  equation  (75)  gives 


0fc+ 1  —  Ok  +  I{0) 


-i  dlog[L(y,  d)\ 
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which  is  the  method  of  scoring.  The  method  of  scoring  is  presented  by 
a  derivation  of  the  approximate  diagonal  information  matrices  for  the 
parameterization  with  u  and  v.  These  results  are  presented  with  different 
notation  than  that  for  amplitude  and  frequency,  which  has  proven  to  be 
useful  for  implementation  in  MATLAB. 

For  the  single  harmonic  case  s(u,  v,  u>)  =  Cu  +  Sr>,  so  that 


ds 

dd  ~  1 

[C  S  ' 

(77) 

where 

ds 

dco 

dC  ds 

s^n+ezv 

(78) 

and 

rdci 

.du>\i 

— t(i )  sin[cj£(i)]  , 

rdSl 
du\ i 

t(i)  cos [uit(i)\ 

(79) 

are  nxl  vectors  of  first  derivatives  obtained  from  equation  (55).  The  infor¬ 
mation  matrix  is  approximately  diagonal 


r  cTc  o  cTf5 1  ho  o' 

I(u,v,u)  =  ^  o  STS  Sr^  ~  ^2  0  1  0  ,  (80) 

"  [oo  &AW. 

which  greatly  simplifies  the  computation  for  each  iteration  in  equation  (76). 
For  multiple  harmonics  with  0  =  (u,  v,  oo), 


ds 

ae 


c  s 


ds  ds 

du>\  ' ' '  dujm  J  ’ 


(81) 


where 


ds  dCk  dSk 

n —  —  a - Uk  +  a - vk 

duk  aook  du>k 


(82) 


and  generalizes  equation  (77).  The  information  matrix  is  still  approximately 
diagonal, 

10  0 

Ti 

/(u,v,u?)«^2  oi  0  )  (83) 

Lo  0  n2 D/12. 


where  D  is  given  in  equation  (27),  which  once  again  leads  to  efficient  com¬ 
putation  of  each  iteration  in  equation  (76).  The  Hessian  for  the  harmonic 


case  is  similar  to  that  of  the  multiple,  except  that  it  contains  only  one  fre¬ 
quency  parameter  instead  of  m  frequency  parameters.  With  6  =  (u,  v,  u), 

can  be  written  as  in  equation  (77),  but,  in  this  case. 


ds  _  dX  _  dC  ds 
dio  duj  duj  ^  du) 


(84) 


and 


SC 

duo 

SS 


Sw 


*? 


n 


=  —hjt(i)  sin {hjUt(i)\  , 
=  hjt.(i)  cos  [hjut(i)\ 


(85) 


are  n  x  m  matrices  of  first  derivatives  obtained  from  equation  (55),  with 
ujk  =  hktjo.  The  information  matrix  is  once  again  approximately  diagonal. 


/( u,  v.w) 


n 

2cr* 


10  0 

0  1  0 

[o  0  J^n2  h2A2  J 


(86) 


and  leads  to  an  efficient  implementation  of  equation  (76)  to  estimate  fun¬ 
damental  frequency. 

As  expected,  the  information  matrices  in  equations  (80),  (83),  and  (86)  lead 
to  the  same  bounds  for  the  frequency  parameters  in  equations  (17),  (28), 
and  (30),  respectively.  One  can  also  get  the  approximate  CRBs  for  u  and  v. 


CRB[uk\ 
CRB[vk ] 


2c? 

n 

2 ^ 

n 


(87) 


for  both  the  multiple  and  harmonically  related  models.  The  CRBs  in  equa¬ 
tion  (87)  are  obtained  exactly  for  the  case  when  the  frequency  parameters 
are  known. 

Because  of  equation  (57),  methods  can  be  developed  for  only  the  MLE  of 
the  frequency  parameters.  For  the  general  multiple  harmonic  case,  this  is 
equivalent  to  finding  cu  to  minimize 

(w)  =  SS( u,  v, u>)  =  eTe  =  ||y  -  Py ||2  =  ||y ||2  -  yTPy  ,  (88) 

or,  alternatively,  to  maximize  the  signal  energy 

E{<jj)  =  yrPy  ,  (89) 
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where 


P  =  PX  =  X(XTX)_1XT  (90) 

is  the  projection  onto  the  subspace  spanned  by  the  columns  of  X(u>).  Note 
that  e,  defined  in  equation  (88),  is  the  residual  error  at  the  optimal  value. 

Newton's  method  for  estimating  fundamental  frequency  is 


fdSS2\  1  OSS 
^  d2u  )  du > 


(91) 


where  SS  =  SS(u)  is  given  by  equation  (88)  for  the  harmonically  related 
case.  To  simplify  the  implementation  of  equation  (91),  I  propose  a  decou¬ 
pled  procedure  that  first  estimates  the  coefficients  assuming  fixed  frequency, 
perhaps  using  an  efficient  approximation  such  as  equation  (68),  and  then 
optimizes  frequency  assuming  fixed  coefficients.  It  turns  out  that  this  pro¬ 
cedure  is  equivalent  to  setting  the  cross-terms  in  the  information  matrix 
involving  (u,  v)  and  ui  to  zero  in  the  method  of  scoring  procedure  for  the 
coupled  harmonic  case.  For  example,  the  method  of  scoring  for  this  ap¬ 
proach  simply  becomes 


where,  as  mentioned  before,  ft  can  be  an  approximate  MLE  of  the  coeffi¬ 
cients  to  conserve  computation.  It  is  important  to  note  that  defining  the 
phase  in  the  middle  of  the  time  window  is  not  just  a  matter  of  convenience 
in  this  case,  but  is  absolutely  necessary  for  the  performance  of  this  decou¬ 
pled  estimation  procedure. 
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6.  Colored  Noise 


Now  I  consider  the  case  where  Z  is  no  longer  assumed  white,  but  instead 
wide-sense  stationary  (WSS)  Gaussian  with  a  known  covariance  of  Kzz(k ) 
and 

a2(u)  =  PZzH=  £  e~i“kKzz(k)  (93) 


is  the  power  spectral  density  or  local  noise  variance.  If  the  noise  vector  z  in 
equation  (8)  is  multivariate  Gaussian  Nn(0 ,  E),  where  Ej j  =  Kzzij  —  j), 
then 

l^  =  %^Te  (94) 


is  the  information  matrix.  Approximate  results  in  the  colored  Gaussian 
noise  case  can  be  obtained  with  the  following  approximation,  which  is 
valid  for  large  enough  n  and  under  mild  assumptions  on  a2  (to) 


arE_1b 


A(to)B*(to)  . 

~duJ  ’ 


27r  J — 7r  <r2(u;) 


where  a  and  b  are  n  x  1  vectors  with  continuous  Fourier  transforms  (i.e., 
equation  (67))  A(to)  and  B( to),  respectively.  With  the  use  of  equation  (95)  for 
the  case  of  estimating  multiple  and  coupled  harmonics,  the  vectors  often 
asymptotically  have  discrete  spectra,  and  the  approximation  can  be  eval¬ 
uated  quickly  by  the  notion  of  delta  functions  or  more  rigorously  by  the 
measure-theoretic  interpretation  of  the  integral.  An  important  point  to  re¬ 
alize  about  this  approximation  is  that,  depending  on  the  nature  of  the  co- 
variance  function,  it  may  require  a  larger  n  to  become  valid  than  is  given  in 
equation  (47). 


An  interesting  result  of  equation  (95)  (or  methods  used  by  Ghogho  and 
Swami  [8]),  is  that,  under  mild  assumptions,  the  CRBs  for  the  colored  Gaus¬ 
sian  noise  case  are  approximately  equal  to  those  in  equations  (28)  and  (30), 


the  local  SNR.  The  dependence  on  local  SNR  can  be  seen  heuristically  by 
noting  that  the  spectral  density  for  the  partial  derivative  of  s  with  any  of 
the  harmonic  signal  model  parameters  will  asymptotically  have  all  its  mass 
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concentrated  at  the  corresponding  frequency,  and,  by  equation  (95),  the  en¬ 
try  in  the  information  matrix  will  depend  only  on  the  local  noise  variance. 
The  MLE  for  colored  Gaussian  noise  is 

P  =  (XrE_1X)_1XrE_1y  ,  (97) 

which  generalizes  equation  (57).  Using  equation  (95),  it  can  be  shown  that 

CTE_1S  ft!  0 

StE~1C  ft!  0  ,  (98) 

and  equation  (97)  separates  into  the  cosine  and  sine  components  and  in  the 
same  way  as  equation  (60).  Further, 

CrE-1C  «  StE-1S  ft!  |diag[lA72(Wl),  1/<72(u;i),  . . . ,  l/<r2(u;m)]  (99) 
and 

cf  E-V  ft  CTky /a2(cok)  , 

S^E-V  ft  S£y /a\uk)  ,  (100) 

so  that 

uk  ft  [2 cr2(u>fc)/n]  [c£y/cr2(u;fc)]  =  (2/n)C^y  , 

vk  ft  2cr2(uik)/n^  [s^y/u2(u;fc)  =  (2/n)S£y  ,  (101) 

which  is  exactly  the  same  approximation  as  equation  (65)  and  the  CFT 
estimates,  equation  (70).  In  a  similar  manner,  the  DFT  estimates  in  equa¬ 
tion  (68)  and  the  estimates  in  equation  (66)  are  approximately  valid  in  the 
colored  noise.  So  asymptotically,  estimating  multiple  harmonic  in  colored 
noise  completely  decouples  into  estimating  m  single  harmonics  in  white 
noise  with  variance  equal  to  the  local  noise  variance  at  that  harmonic. 
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7.  Coupled  Harmonic  Model  for  Battlefield  Acoustic  Targets 


Target  identification  using  battlefield  acoustic  sensor  arrays  is  an  impor¬ 
tant  problem  for  the  Army  [9].  Acoustic  signatures  of  targets  of  interest, 
such  as  tanks  and  trucks,  exhibit  time-varying  patterns  of  harmonic  am¬ 
plitudes  that  facilitate  target  ID.  The  spectrogram  from  one  of  an  array 
of  seven  acoustic  sensors  taken  by  the  Army  Research  Laboratory  (ARL) 
during  recent  field  tests  at  Aberdeen  Proving  Ground  (APG),  Maryland, 
is  shown  in  figure  1.  The  sampling  rate  is  T  =  2000  samples /s  and  FFTs 
were  taken  with  n  =  2048  samples.  The  target  is  an  M60  tank  that  is  ap¬ 
proaching  the  sensor  until  the  time  of  the  closest  point  of  approach  (CPA), 
at  approximately  t  =  40  s,  and  retreating  afterwards.  The  range  at  CPA  is 
approximately  50  m. 

The  harmonic  structure  and  the  time-varying  nature  of  the  signal  are  very 
apparent.  The  strongest  line  at  approximately  90  Hz  is  the  sixth  harmonic, 
so  that  the  fundamental  frequency  is  approximately  15  Hz.  A  coupled  har¬ 
monic  signal  model  accurately  models  a  large  portion  of  a  vehicle's  acoustic 
signature,  especially  that  coming  from  the  engine.  For  the  current  ARL  tar¬ 
get  ID  system  [10],  the  second  through  twelfth  harmonics  are  used  so  that 


Figure  1.  Battlefield 
acoustic  data  with 
time-varying  coupled 
harmonic  signal. 


Time  (s) 
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h  =  (2, 3, 4. 5, 6, 7, 8, 9. 10, 11, 12)r  and  m  =  11  in  equation  (6).  The  addi¬ 
tive  noise  process  Z  in  equation  (7)  consists  of  the  cumulative  effect  of  am¬ 
bient  noise  such  as  wind,  broadband  and  narrowband  nontarget  sources, 
and  a  portion  of  the  target's  acoustic  signature  not  completely  described 
by  equation  (6),  such  as  the  additional  line  at  approximately  50  Hz  in  fig¬ 
ure  1,  which  is  coming  from  the  track  slap.  An  example  of  data  from  a  short 
time  window  (approximately  0.25  s)  is  shown  in  figure  2,  which  shows  the 
first  512  samples  from  the  same  data  as  in  figure  1.  For  battlefield  acoustic 
scenarios,  the  fundamental  frequency  can  be  assumed  constant  for  short 
time  windows  (on  the  order  of  1  s),  and  the  stationary  methods  presented 
can  be  applied. 


7.1  Fundamental  Frequency  Estimators  and  Trackers 

For  each  fundamental  frequency  /,  the  residual  error  of  the  coupled  har¬ 
monic  signal  model  can  be  efficiently  and  simply  estimated  by  the  FFT 
of  y,  with  fk  =  hkLo/( 2tt)  in  equation  (68).  In  practice,  the  search  for  the 
fundamental  frequency  is  limited  over  a  small  range  of  values,  e.g.,  5  to 
20  Hz  in  increments  of  0.1  Hz.  A  plot  of  the  residual  error  power  SS/n  in 
equation  (88)  for  the  data  in  figure  2  (the  true  fundamental  frequency  is  ap¬ 
proximately  /  =  15.5  Hz)  is  shown  in  figure  3.  A  series  of  clear  regions  of 
minimum  power  shows,  in  this  case,  that  the  fundamental  frequency  needs 
to  be  known  initially  to  within  a  certain  accuracy  (roughly  5  Hz)  to  guaran¬ 
tee  that  this  method  finds  the  correct  local  minimum.  Once  a  fundamental 


Figure  2.  Acoustic  data 
with  coupled  harmonic 
signal. 
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Figure  3.  Residual  error 
power  for  FFT 
fundamental  frequency 
estimation. 


Figure  4.  Residual  error 
power  for  MLE 
fundamental  frequency 
estimation. 


frequency  is  estimated,  the  estimate  of  amplitude  can  be  readily  found,  al¬ 
though  it  is  important  to  use  the  corrected  estimate  of  the  phase  in  equation 
(69).  The  fundamental  frequency  can  be  estimated  more  accurately  with  the 
exact  MLE.  A  plot  of  the  residual  error  power  for  the  MLE  of  the  data  in 
figure  2  (assuming  white  Gaussian  noise)  is  shown  in  figure  4.  The  plot 
is  similar  to  figure  3,  except  that  the  maximum  values  are  much  sharper, 
which  suggests  that  the  MLE  estimate  will  be  much  more  accurate. 
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Figures  3  and  4  suggest  that  a  simple  but  effective  fundamental  frequency 
tracker  can  be  constructed,  assuming  an  accurate  initial  value  has  been  ob¬ 
tained.  Subsequent  estimates  can  be  obtained  if  one  minimizes  equation  (3) 
over  a  range  of  fundamental  frequencies  equal  to  the  previous  estimate 
plus  a  reduced  number  of  offsets  reflecting  the  maximum  expected  change 
over  the  time  window.  For  the  data  in  equation  (1),  accurate  estimates  of 
the  fundamental  frequency  were  efficiently  obtained  using  FFTs  with  this 
approach,  an  initial  value  of  15.5  Hz,  offsets  of  —0.25  to  0.25  Hz  in  incre¬ 
ments  of  0.025  Hz,  and  time  windows  of  length  n  =  2048  with  75  percent 
overlap.  The  iterative  procedures  using  the  gradient  descent  methods  de¬ 
scribed  previously  can  be  modified  slightly  to  track  fundamental  frequen¬ 
cies  by  an  adaptive  filter.  For  the  data  in  equation  (1),  this  approach  tracked 
the  fundamental  frequency  accurately  for  most  of  the  run,  but  did  lose  the 
signal  for  the  last  10  s  or  so,  where  the  signature  changed  significantly. 

Harmonic  Amplitudes  and  Phases  for  Target  Identification 

Using  the  estimates  of  fundamental  frequency  obtained  from  the  tracker 
described  above,  I  obtained  estimates  of  amplitude  and  phase  for  the  data 
in  figure  1  using  equation  (60).  Figure  5  shows  the  estimated  amplitudes  of 
the  third  and  sixth  harmonics;  despite  the  variation  in  power  as  a  function 
of  range  due  to  acoustic  propagation  effects,  these  amplitudes  have  a  fairly 
constant  difference.  This  suggests  that  the  ratio  of  the  third  and  sixth  har¬ 
monics  (the  strongest  line)  is  a  more  robust  feature  for  target  ID.  For  this 
reason,  normalizing  the  amplitudes  by  the  strongest  line  is  part  of  ARL's 
current  target  ID  system. 

A  phase  coupling  also  exists  between  the  third  and  sixth  harmonics,  as 
shown  in  figure  6,  which  is  a  plot  of  the  difference 

dse  =  mod(2 03  -  06,  2w)  (102) 

versus  time.  The  value  of  g?36  is  relatively  constant  over  short  time  win¬ 
dows.  However,  the  difference  is  not  zero,  as  predicted  in  some  models 
(e.g.,  see  Swami  and  Ghogho  [7]),  and  it  does  appear  to  vary  slightly  over 
the  entire  time  period.  This  apparent  anomaly  is  perhaps  due  to  propa¬ 
gation  effects  and/or  variations  in  the  target's  aspect  or  operation  char¬ 
acteristics.  The  extent  to  which  differences  such  as  in  equation  (102)  are 
useful  in  target  identification  remains  an  open  issue,  as  does  the  accuracy 
with  which  these  differences  can  be  estimated.  One  reason  the  harmonic 
numbers  3  and  6  were  chosen  to  show  coupling  is  that  the  simple  2  to  1 
ratio  between  them  minimizes  the  variance  of  the  corresponding  differ¬ 
ence.  With  harmonics  5  and  7,  for  example,  the  corresponding  difference 
is  (I57  =  mod ( 7</>5  —  507, 27r)  which,  with  all  else  being  equal,  would  have  a 


Figure  5.  Harmonic 
amplitude  estimates. 


Figure  6.  Phase  coupling 
between  third  and  sixth 
harmonics. 


variance  18.5  times  that  of  equation  (102).  So,  for  targets  where  5  and  7  are 
dominant  harmonic  lines,  the  harmonic  coupling  will  be  more  difficult  to 
detect  and  track  and  may  not  be  a  stable  feature. 

The  importance  of  using  the  corrected  phase  in  equation  (69)  can  be  demon¬ 
strated  in  the  calculation  of  equation  (102).  Without  the  correction,  empiri¬ 
cal  results  showed  that  half  the  time  the  calculated  differences  were  about 
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180°  off  (i.ev  the  worst  possible  error),  and  were  nearly  correct  the  other 
half.  The  reason  for  this,  assuming  uniform  distribution  within  the  Fourier 
frequency  bins,  is  that  half  the  time  the  errors  in  the  phase  estimates  are  in 
opposite  directions  and  add  together  destructively  in  equation  (102),  and 
half  the  time  they  cancel  each  other  out. 

For  large  n,  the  variance  of  the  MLEs  of  amplitude,  phase,  and  fundamen¬ 
tal  frequency  are  approximately  equal  to  the  CRB.  Looking  at  equation  (30), 
the  variance  of  fundamental  frequency  goes  down  with  the  presence  of 
high  SNR  tones  with  large  harmonic  numbers.  For  battlefield  acoustics,  this 
is  not  always  the  case,  because  increased  propagation  loss  occurs  at  higher 
frequencies.  Also,  the  fundamental  frequency  can  be  estimated  fairly  ac¬ 
curately  with  small  data  sizes,  but  equation  (21)  shows  that  exploiting  the 
harmonic  phase  coupling  for  target  ID  may  require  more  time  samples. 
However,  because  of  the  nonstationarity  of  the  fundamental  frequency,  it 
is  not  always  possible  to  increase  n  arbitrarily,  and  the  accuracy  of  phase 
estimates  is  limited  with  stationary  approaches.  So,  obtaining  more  accu¬ 
rate  phase  estimates  with  time-varying  fundamental  frequency  is  a  subject 
for  further  research  and  presents  an  interesting  problem  in  nonstationary 
signal  processing. 


8.  Conclusions 


Efficient  methods  for  the  maximum  likelihood  estimation  of  parameters  for 
multiple  and  coupled  harmonics  were  presented  that  exploit  the  diagonal- 
ity  of  the  information  matrices  when  the  phase  is  defined  in  the  middle  of 
the  time  window.  With  the  generalization  to  local  noise  variances  and  local 
SNRs,  the  CRBs  for  the  colored  noise  case  are  asymptotically  equivalent  to 
the  white  noise  case,  perhaps  requiring  a  larger  number  of  samples  n.  With 
the  same  caveat,  the  approximate  MLE  methods  for  the  white  noise  case  are 
also  applicable  to  the  colored  noise  case.  The  MLE  methods  were  applied 
to  real  battlefield  acoustic  data  and  the  problem  of  target  identification.  The 
difference  between  the  power  levels  of  strong  harmonics  was  shown  to  be 
fairly  constant  over  the  entire  run,  which  demonstrates  the  potential  ro¬ 
bustness  of  these  features.  Harmonic  phase  coupling  was  also  shown  to  be 
present  in  these  acoustic  signatures,  although  challenges  remain  to  fully 
exploit  these  relationships  for  target  ID.  For  example,  with  all  else  being 
equal,  frequencies  with  simple  harmonic  relationships,  such  as  2  to  1  or  3 
to  2,  are  preferred  to  produce  stable  features.  Also,  naive  use  of  FFT  ap¬ 
proaches  to  estimate  phase  leads  to  erratic  estimates.  The  CRBs  presented 
provide  insights  into  the  limitations  of  the  accuracy  of  phase-coupling  esti¬ 
mates  and  suggests  the  value  of  developing  nonstationary  methods. 
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